High-precision determination of the light-quark masses from realistic lattice QCD 
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Three-flavor lattice QCD simulations and two-loop perturbation theory are used to make the 
most precise determination to date of the strange-, up-, and down-quark masses, m s , m u , and md, 
respectively. Perturbative matching is required in order to connect the lattice-regularized bare-quark 
masses to the masses as defined in the MS scheme, and this is done here for the first time at next-to- 
next-to leading (or two-loop) order. The bare-quark masses required as input come from simulations 
by the MILC collaboration of a highly-efficient formalism (using so-called "staggered" quarks), with 
three flavors of light quarks in the Dirac sea; these simulations were previously analyzed in a joint 
study by the HPQCD and MILC collaborations, using degenerate u and d quarks, with masses as 
low as m s /8, and two values of the lattice spacing, with chiral extrapolation/interpolation to the 
physical masses. With the new perturbation theory presented here, the resulting MS masses are 
mf 5 (2 GeV) = 87(0) (4) (4) (0) MeV, and m m (2 GeV) = 3.2(0) (2) (2) (0) MeV, where m = \(ra u + 
ma) is the average of the u and d masses. The respective uncertainties are from statistics, simulation 
systematics, perturbation theory, and electromagnetic/isospin effects. The perturbative errors are 
about a factor of two smaller than in an earlier study using only one-loop perturbation theory. Using 
a recent determination of the ratio m u /md = 0.43(0) (1)(0) (8) due to the MILC collaboration, these 
results also imply m^ 5 {2 GeV) = 1.9(0) (1)(1) (2) MeV and mf S (2 GeV) = 4.4(0) (2) (2) (2) MeV. A 
technique for estimating the next order in the perturbative expansion is also presented, which uses 
input from simulations at more than one lattice spacing; this method is used here in the estimate 
of the systematic uncertainties. 



I. INTRODUCTION 

The strong sector of the standard model contains a 
number of inputs that are a priori unknown and must 
be determined from experiment. Knowledge of these fun- 
damental parameters, the quark masses and the strong 
coupling, also requires precise theoretical input, because 
of confinement: quarks and gluons cannot be observed as 
isolated particles, hence one can only extract their prop- 
erties by solving QCD for observable quantities such as 
hadron masses, as functions of the quark masses and the 
strong coupling. Precise values of the quark masses in 
particular are valuable in many phenomenological appli- 
cations, such as in placing constraints on new physics 
beyond the standard model. 

Recent breakthroughs in lattice QCD are having a sig- 
nificant impact on the determination of many hadronic 
quantities. These advances are due to two related devel- 
opments: the use of perturbation theory in the design 
of so-called improved-lattice discretizations p], Q , which 
significantly reduces cutoff effects; and numerical sim- 
ulations of a highly-efficient formalism (using so-called 
"staggered" quarks), where the correct number of light - 
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quark flavors can be included in the Dirac sea, and at 
sufficiently small quark masses Q , enablingaccurate chi- 
ral extrapolations to the physical region 0, [j| . These 
developments have eliminated the large and frequently 
uncontrolled systematic uncertainties inherent in most 
other lattice QCD studies, where the effects of sea quarks 
are either completely neglected (working in the so-called 
"quenched" approximation), or where the simulations are 
done with the wrong number of sea quarks (and at very 
heavy masses) : comparison of the predictions of such un- 
realistic theories with experiment reveals typical incon- 
sistencies of 10-20%, which precludes their use in many 
important phenomenological applications. 

Recent work by our group, the High Precision QCD 
(HPQCD) collaboration, together with the Fcrmilab and 
MILC collaborations, has instead utilized unquenched 
simulations with up, down and strange quarks in the 
Dirac sea. These simulations are made computationally 
feasible by using the 0(a 2 )-improved staggered-quark ac- 
tion (a is the lattice spacing), with the one-loop 0(a 2 )- 
improved gluon action (hereafter collectively referred to 
as the "asqtad" action @,|3|)- We have shown that these 
three-flavor simulations give results for a wide variety of 
observables that agree with experiment to within sys- 
tematic errors of 3% or less p. This precision, and the 
unified description of hadronic physics in many different 
systems, is a striking consequence of having a realistic 
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description of the Dirac sea. 

The application of this program to many other quan- 
tities of interest requires additional use of perturbation 
theory, in order to properly match lattice discretiza- 
tions of the relevant operators and couplings onto con- 
tinuum QCD. These perturbative-matching calculations 
must generally be carried out at next-to-next-to-leading 
order ("NNLO" order, which is generally equivalent to 
two-loop Feynman diagrams), if one is to obtain results 
of a few-percent precision 0, 01 ■ The success of this ap- 
proach has been demonstrated by a recent HPQCD de- 
termination of the strong coupling at NNLO order, with 
the result a^jg(M z ) — 0.1170(12) jgj, the most accurate 
of any method. 

In this paper we apply our perturbation theory pro- 
gram to make the first-ever NNLO determination of 
the light-quark masses (the computationally-simpler case 
of the additive zero-point renormalization for Wilson 
fermions was previously computed at two-loop order in 
Ref. 0| , and for static quarks in Ref. ^3). The quark 
masses are not physically measurable, and as such are 
only well defined in certain renormalization schemes, 
such as the MS mass m MS (^), evaluated at some rele- 
vant scale fi. We have done the necessary perturbative 
matching calculations to connect the lattice-regularized 
bare mass mo (a) for the "asqtad" action, as a function 
of the lattice spacing, to m (/x). 

The lattice spacings and bare quark masses are re- 
quired as input, and these have been determined in earlier 
studies. The lattice spacings we use were determined ul- 
timately from the T' — T mass difference p^. which is 
insensitive to the quark masses; the results agree within 
systematic errors of 1.5-3% with the lattice spacing ex- 
tracted from a wide variety of other physical quantities 
; this eliminates an uncontrolled systematic error in 
studies that are done without the correct number of sea 
quarks. 

The bare quark masses have previously been deter- 
mined in a joint HPQCD and MILC collaboration study 
[}M- using partially-quenched chiral-perturbation theory 
[3, |5| to make precise extrapolations to the physical re- 
gion; the simulations used equal dynamical u/d quark 
masses as small as m s /8. These bare masses were used 
in Ref. 0] to estimate the MS masses using one-loop 
perturbation theory [l4j |; the dominant systematic er- 
ror in that determination came from unknown second 
and higher orders in the perturbative matching. Sig- 
nificant progress is made here due to our computation 
of the second-order perturbative matching coefficient. 
We thereby reduce the systematic error from pertur- 
bation theory by about a factor of two compared to 
the earlier result of Ref. The remaining perturba- 

tive error is the same size as the current lattice system- 
atic error, which is largely due to the chiral extrapola- 
tion/interpolation. 

The staggered-quark formalism had previously been af- 
flicted with several poorly-understood problems, which 
have been tamed with an aggressive program of pertur- 



bative improvement, leading to the "asqtad" action that 
we are using here. With an unimproved-staggered ac- 
tion large discretization errors appear, although they are 
formally only 0{a 2 ) or higher. Many of the renormaliza- 
tion factors required to match lattice operators onto con- 
tinuum quantities are also large and poorly convergent 
in perturbation theory for unimproved-staggered quarks; 
this is true, for example, of the mass renormalization 
that is needed here. It turns out that both problems 
have the same source, a particular form of discretiza- 
tion error dubbed "taste violation," and both are ame- 
liorated by use of the improved-staggered formalism 0. 
The perturbation theory then shows small renormaliza- 
tions |^I^ JT3 . ll5llT^ |. and discretization errors are much 
reduced Qj- Taste violation is strongly probed in cer- 
tain quantities, notably the would-be Goldstone meson 
masses which play an important role in our analysis, 
and these effects are taken into account in the staggered 
chiral-perturbation theory £| . 

A potentially more fundamental concern about stag- 
gered fermions relates to the need to take the fourth 
root of the quark determinant, in order to convert the 
four-fold duplication of "tastes" into one quark flavor. 
One might imagine that the fourth root would intro- 
duce nonlocalities (and would induce violations of uni- 
tarity), which would prevent decoupling of the ultravio- 
let modes of the theory in the continuum limit. However 
a great deal of theoretical evidence has been amassed 
which demonstrates that the properties of the staggered 
theory, with the fourth root, are equivalent to a one- 
flavor theory, up to the expected discretization errors [|| ; 
these are again due to short-distance taste-changing in- 
teractions, which are mediated by high-momentum glu- 
ons [2| . The locality of the free-field staggered theory is 
trivial, and is made manifest in the "naive" basis used 
in Ref. Q. Non- localities do not arise in perturbation 
theory, since the staggered-quark matrix is diagonal in 
the taste basis, up to those short-distance (and calcula- 
ble) corrections. Most important, it has been shown that 
perturbative improvement of staggered actions correlates 
exceedingly well with nonperturbatively-measured prop- 
erties of the staggered-fermion matrix, providing clear 
support for the correctness of the fourth-root procedure. 
This includes the measured pattern of low-lying eigenval- 
ues of the staggered matrix [l7j, and the measured pat- 
tern of taste-violating mass differences in the non-chiral 
pions [is| . 

The rest of this paper is organized as follows. Sec- 
tion [n] details the computation of the two- loop matching 
factor, using the pole mass as an intermediate match- 
ing quantity. In Sect. IIIII we use the bare-quark masses 
from the MILC "asqtad" simulations to extract the MS 
masses, including an analysis of the systematic uncer- 
tainties. In that connection, we also describe how our 
perturbation theory results can be used to estimate the 
third order in the perturbative expansion, using input 
from simulations at more than one lattice spacing. Sec- 
tion presents some conclusions, including a compari- 
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son of our results with other recent determinations of the 
strange-quark mass. 

II. PERTURBATION THEORY 
A. Lattice to MS matching 

We do two-loop perturbative matching to connect the 
cutoff-dependent lattice bare-quark masses mo (a) to the 
MS masses m MS (/i) at a given scale /i. We define the 
perturbative matching factor Z m ([ia,moa) according to 

m °(n) = Z m (fxa,m a), (1) 

a 

where we make explicit the fact that the simulation input 
provides the bare masses in lattice units. We compute 
Z m in two stages, using the pole mass M as an inter- 
mediate matching quantity. We also use our previous 
determination pi Il6j of the relation between the lattice 
bare coupling aq a t and the renormalized coupling ay(g), 
defined by the static potential, to reorganize both sides 
of the matching equation into series in terms of ay {q* ) , 
at an appropriately-determined scale. 

In the following subsections we consider in turn the 
pole mass in the MS- and lattice-regularization schemes. 
We also provide some details on the consistent evaluation 
of the two-loop on-shcll condition in the lattice scheme, 
and the various checks that we have applied to our eval- 
uation of the two-loop self-energy. We then quote our 
results for the matching factor of Eq. QJ, and for the 
relevant matching scale q*. 



and where the contribution zfh from an internal-quark 
loop with the same flavor as the valence quark is split off 
from the contribution zfl of ni internal-quark loops with 
different flavor (these are taken here to be degenerate 
in mass, though this is trivially generalized). The total 
number of flavors is rif — ng + 1. The individual functions 
are given by 

Zp = — 1 — f^/xM, 

Z FF = 128 ~ T^ 2 _ + 3 ^2 log 2 + 32^M + w^M) 

— 1111 _i_ 1 r _i_ 3 /■ V l^rrO 185tf 11«2 
Z FA - -"384" + 2^2 + s<=>3 ~ 2& l0 8 2 ~ "9T^ M ~ 32V*f' 

zfl = H + \(,2 - 2A(r sca ) + ^I^m + §^?m> 

zfh = i + |Ca - 2A(1) + if e^M + \^M, (4) 

where = \og{^ 2 /M 2 ), 

m sca . . 

''sea = , (p) 

m V alcnce 

and where m sca and ra va i e nce are the sea- and valence- 
quark masses, respectively. The function A(r) gives the 
dependence of the renormalization factors zfl and zfh 
on the quark mass in an internal-quark loop (sea and 
valence, respectively). An exact integral expression for 
A(r) can be found in |2Cj : particular limits are 

A(r«l) = |C 2 r + 0(r 2 ), 



A(r > 1) = |log 2 r + if log r 

+ 1C2 + M+°( r ~ 2 log r )- ( 6 ) 



B. Pole mass in the MS scheme 

We begin by recalling the relation between the MS 
mass and the pole mass which is known 

through three loops [2ll |2^ . We require it to second 
order, a result that was first obtained in Ref. [2(j (ex- 
pressions for the relation at arbitrary fj, are conveniently 
given in Ref. plf) 



,MS 



(ji) = M 



+ Z2 



1 + Zl 



(-) 



\MJ 7T ; 



(m) 



(2) 



where the one- and two-loop coefficient functions 
z%(fj,/M) and Z2(fi/M) are reduced to a set of terms 
with different color structures [in the following Cf — 
(iV 2 — l)/(2iV c ), Ca = N c , and T — 1/2] 



z\ = C F z F , 

Zl = CpZFF 



CfGazfa + C F Tri£ZFL + C f Tz F h, 



(3) 



C. Pole mass in the lattice scheme 

The relation between the pole mass M and the bare 
quark mass mo (a) in the lattice-regularized theory has 
the form 



M = m 1 + ai a t (.An log m a + A w ) 

+ af at (A 2 2 log 2 m a + A 2 \ log m a + A 2 < 



2(1 ) 



(7) 



where ai at (a) is the bare lattice-regularized coupling. 
The coefficients of the logarithms are determined by the 
anomalous dimension in mo (a), which in turn can be de- 
termined from the known anomalous dimension of the 
MS mass, as described in Ref. |2j|. This also requires the 
connection between ai a t and ajjg, which can be found 
from our evaluation of the relation between ct\ a t and 
a\r(q), which we obtained at NNLO order in Refs. 
(see also Ref. [24]); here we only require the NLO relation 



aiat = av{q) [1 - v\{q)av{q) 



(8) 
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where 

Vl (q) = ^\n(—) +v w , /3 = ll-^ / , (9) 
47r \ aq J 3 

and where the constant vio for the "asqtad" action is 

uio = 3.57123(17) - 1.196(53) x 10~ 4 n f . (10) 

We then obtain (SU(3) color is used throughout) 

An = 2/tt, 

A 22 = 0.7599 - 0.03377 n/, 

A21 = -4.2164 + 0.09429 n f - 0.6366 A w . (11) 

The terms A\q and A20 require the computation of the 
one- and two-loop lattice self-energies, respectively. We 
have previously determined the first-order term for the 
"asqtad" action (HQ 

(12) 



A w = 0.5432(1), 



neglecting corrections of O((am ) 2 ). What is new here 
is the evaluation of ^201 which is detailed in the next 
subsection. 

As with the quark-loop terms in the MS connection to 
the pole mass, A 2 o depends on the number of quark fla- 
vors in the sea, and on the loop-quark masses. However, 
as we demonstrate explicitly in the next subsection, the 
mass dependence in A20 cancels precisely in the matching 
to the continuum relation Eq. J2J), in the limit where the 
sea and valence quark masses are all much less than the 
lattice ultraviolet cutoff n/a, and the MS scale fx. The 
cancellation holds for an arbitrary ratio m sea /m va i ence 
[compare with the mass-dependent corrections in the con- 
tinuum quantities zfl and Zfh in Eqs. © and fllf]. This 
follows from the fact that, in the limit where the energy 
scales are large, the net-renormalization factor connect- 
ing mo (a) to m MS (p) probes internal loops at scales large 
compared to the quark masses; hence the dependence on 
quark masses of the intermediate renormalization factors 
connecting to the pole mass M are infrared effects that 
are identical in the lattice- and MS-regularized theories. 



D. Evaluation of the two-loop lattice self-energy 

The renormalizcd pole mass M is identified from the 
pole in the dressed-quark propagator G(p, mo), 



G(p, m ) 1 = if+ m + £tot(p), 



(13) 



where £ to t is minus the usual 1PI two-point function. 
The two-loop diagrams contributing to S to t in the lat- 
tice theory are shown in Fig. ^ The lattice-dispersion 
relation implied by p will be kept as a general function of 
the lattice momentum p; for improved-staggered quarks, 



for instance, ap^ = sin(ap M ) 
spinor decomposition: 



i sin 3 (ap p ). We make the 



£ to t(p) = + + m o) £ 2 (p), 



(14) 



where Si and S2 are both implicit functions of mo, and 
both are Dirac scalars. At zero three-momentum the 
renormalized on-shell condition is given by 



P(Pt) 



-ipt = m 



Sib) 



1 + S a (p)' 

the solution of which determines the pole mass 
Pt = %M. 

The perturbative expansion of M is denoted by 

M = M (0) 



(15) 



(16) 



(17) 



with a corresponding notation for other quantities (in- 
cluding, e.g., p[ 0) = iAf(°\ where P(iM^) = m ). At 
first order one has 



M( 2 > = 4 0) £fV°>), 



where 



Z. 



(or 



dP{ix) 



dx 



(18) 



(19) 



=M(°) 



is the tree-level wave function residue. 

An evaluation of the on-shell condition to second order 
requires a consistent expansion of the right-hand side of 
Eq. I|15jl. Part of the 0(g$) term arises from the one- 
loop piece of Si(p t ), when it is evaluated at the one- 
loop-corrected on-shell energy 



^\p t )^^\ lM ^) + l9 lM^ 



d^f\pt) 



dp t 



(20) 



Hence the two-loop contribution to the pole mass is given 
by 



M( 4 ) = Z ( , 0) I m< 4 > 



2 d 2 P{ix) 



dx 2 



(21) 



where 

mW = E[ 4) (iM(°') 
with 



EfliMf ')^'^' '), (22) 



,(2) 



(*M(°>) = -S 2 2) (iM(°)) +iZ. 



, (0) asfV 



dp t 



tM(») 

(23) 



Note that the second-term in Eq. Ij21|l is a correction of 
O((amo) 3 ) in the continuum limit of the "asqtad" action. 

The derivative £j Q j! of the one-loop self-energy (differ- 
entiation is implicitly defined in Eq. [12: il) as with respect 
to P{pt)) is also the one-loop part of the wave function 
renormalization (up to finite-lattice discretization correc- 
tions), and is infrared divergent. This infrared divergence 
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FIG. 1: The two-loop diagrams that contribute to the pole-mass renormalization in the lattice scheme. The numbering is 
consistent with Ref. :1JJ] . The filled circles represent five one- loop sub-diagrams which dress the gluon propagator at that order 
(these contain gluon and ghost loops, and the measure term), in addition to the diagrams with internal-quark loops, which 
are shown explicitly. The quark-loop diagrams are understood to be summed over all flavors. The crosses represent vertices 
generated by the perturbative expansion of tadpole and other renormalization factors in the gluon and quark actions. 



precisely cancels against an infrared divergence in the 
two- loop nested-rainbow diagrams (the sum of D21, D22, 
and CT4 in Fig. , which parallels the cancellation of 
infrared divergences in the continuum self-energy at this 
order. In this connection, we note that the continuum- 
pole mass was first shown to be infrared finite at two 
loops in Ref. i an all-orders proof of its finiteness has 
been established only fairly recently, by Kronfeld plj . 

Figure[5]illustrates the infrared cancellation in Eq. i|22|) 
in diagrammatic form. We numerically evaluate the two- 
loop integral for the left-hand diagram in the top line of 
Fig. [21 including, in its integrand, a term with the prod- 
uct of the independent integrands for the one-loop self- 
energy, and its derivative. This grouping is IR finite, 
and does not require any infrared regulator. We obtain 
a stringent check of this result by noting that this com- 



bination generates a leading-logarithmic contribution to 
the anomalous dimension of the mass which goes like 
log 2 (amo), and whose coefficient is identical to that of 
the infrared-subtracted combination in the continuum, 
which can easily be read off in Feynman gauge from the 
MS results in \]J%. 

All the diagrams for the two-loop self-energy, Fig. ^ 
were generated and evaluated independently by two of us. 
The Feynman rules for the highly-improved actions are 
exceedingly complicated, and were generated automati- 
cally using computer-algebra based codes 0, Q • More- 
over, one of us has produced an algorithm which auto- 
matically generates the Feynman diagrams themselves 
[j| . We can then readily evaluate the lattice diagrams for 
a wide variety of gluon and quark actions. The two-loop 
integrals are evaluated numerically, using the adaptive 




= IR finite, 



where 



FIG. 2: Schematic representation of an IR subtraction for the two-loop pole mass; appropriate traces of the self-energy with 
an energy projector are implicit. 



Monte-Carlo method VEGAS 26]. A powerful additional 
cross-check was provided by an explicit verification that 
our results are gauge independent, which we established 
numerically for two different bare-quark masses in three 
covariant gauges: Feynman, Landau and Yennie. 

Knowing the coefficients of the logarithmic terms in 
Eq. (J7J also provides a valuable cross-check on our re- 
sults, and increases the accuracy of our numerical de- 
termination of the remaining term A20 ■ We compute the 
two-loop pole mass as a function of mod, and subtract the 
known logarithms, in order to isolate which must be 




FIG. 3: The n f = part of the two-loop pole mass M (4) 
for the "asqtad" action (in units of a?f at ), for varying quark 
mass moa (upper panel). The squares and crosses distinguish 
two independent sets of calculations, each one by a different 
author. The lower panel shows the results after subtracting 
the known logarithms in M^ A \ which yields the term A20 in 
in Eq. 0. 



finite as niga — > 0. Figure [3] shows the rif = part of the 
results, corresponding to the diagrams without internal- 
fcrmion loops, over a wide range of bare masses, and 
clearly shows the expected limiting behavior. 

Additional checks are provided for diagrams which 
have a leading log 2 (amo) term, which arises from the 
infrared limit of both the outer- and inner-loop integrals, 
and which is therefore an infrared quantity, independent 
of the ultraviolet regulator; the coefficients of these dou- 
ble logarithms in the individual diagrams are available in 
Feynman gauge from the original two-loop MS calcula- 
tion of Tarrach and we have verified that these are 
reproduced in our lattice calculation. 

A further stringent check of our evaluation of the di- 
agrams with fermion loops (diagrams 12, 13, 19 and 20 
in Fig. ^) is achieved by computing over a range of sea 
quark masses. As discussed in Sect. Ill CI the mass de- 
pendence in the intermediate renormalization from the 
bare mass to the pole mass M should cancel against the 
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FIG. 4: Comparison of the loop-quark mass dependence of 
pole-mass renormalization factors, on the continuum side of 
the matching, and on the lattice side. The lattice calculation 
used am V aience = 0.001. The continuum function was taken 
from Ref. j^JJ]; the solid lines are from analytic expressions 
for the associated limits, cf. Eq. (|SJ 
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FIG. 5: BLM scales aq* for the lattice-bare mass to MS mass 
matching factor Z m (^a), as functions of /*a. The solid line 
shows the first-order scale q*, while the dashed and dotted 
lines show the second-order scales gj+ an d q%- , respectively, in 
regions where they are real. The first-order scale is the correct 
choice for fia < 0.61 and fia > 3.74, while the appropriate 
second-order scale applies to the region in between. 



renormalization from M to the MS mass. We define 

4 

A 20 (r S ca) = ^20(0) + =-a Ajat (»•«»), ( 24 ) 

and compare with the analogous continuum function 
A(r sca ), cf. Eqs. We plot our results in Fig. |U 

over a very wide range in r sca , which shows the expected 
agreement. 

Our result for the matching term A 20 , in the limit of 
vanishing sea-quark mass, is 

A 2Q = 6.092(20) - 0.1484(5)^ - 0.0328(3), (25) 

where the last term corresponds to an internal-quark loop 
containing the valence-quark flavor. The uncertainties 
arise from the numerical evaluation of the loop integrals. 

The relatively large value of A 2 q may be a symptom 
of the renormalon ambiguity [23| in the pole mass M; 
there is similarly a large second-order term in the con- 
nection between M and the MS mass, Eq. J2J). However, 
as seen in the next subsection, these large corrections 
in the intermediate matchings to the pole mass mostly 
cancel in the final matching of the bare mass to the MS 
mass, Eq. for which perturbation theory should be 
(and is) reliable. 



E. Results for Z m and its BLM scale 

To complete our determination of the matching factor 
Z m in Eq. we substitute the expression Eq. for 
the pole mass in the lattice scheme into the equivalent 



MS expression in Eq. ©, and expand the logarithms in 
powers of the coupling. We also reorganize the couplings 
to the ay scheme at some scale q*, which we determine 
below according to the BLM procedure j2^. This leaves 
an expression with logarithms only of /ia and aq*. The 
logarithms of mo drop out of Z m , as expected, since these 
are infrared effects that are identical in the intermediate 
lattice and continuum matchings. [We note that the co- 
efficients of the logarithms in Eq. (JTJ), corresponding to 
the lattice anomalous dimension, could instead be deter- 
mined from the requirement that the logarithms of mo 
drop out of the final expression for Z m .\ 
Z m then takes the form 

Z m (na,m, a) = 1 + Z m .i(^a)a v (q*(na)) 

(Ha)a v + O (a 3 , (m a) 2 ) , (26) 

where the first-order term was derived previously in 
Ref. [l| 

^ ro ,i(H = 0.1188(1) - - log(na). (27) 

The new information presented here is the expression for 
Z m ^{^a), which takes the form 

Z m .2(^ a ) = z 22 log 2 (/ia) + Z 2 i log(/ia) + Z 2Q , (28) 
wherre 

Z 22 = 0.75 99 - 0.03 38 n/, (29) 
Z 21 = -0.4049- 0.0281 n f 

+ (-1.1145 + 0.0675 n/)log(aa*), (30) 
Z 2Q = 2.086(20) - 0.0144(5)71/ 

+ (0.2080 - 0.0126 n/)log(ag*). (31) 

To find the optimal choice for aq* in the BLM scheme, 
as a function of /ia, we analytically evaluated the average- 
momentum scales in the one-loop MS self-energy diagram 
(with appropriate renormalizations psf). and made nu- 
merical evaluations of the average scales in the lattice 
self-energy. We find that the second-order scale-setting 
procedure of Ref. [2^| is needed over a wide range of scales 
in the region of /ia = 1. 

At leading order the scale q* is determined by P, |2^| 



log( 0l 2 ) 



/A/(<7)log(g 2 ) _ (/(g)log(q 2 )} 



JdUf(q) 
= «log(<7 2 ))) 



(/(<?)) 



(32) 



where f(q) is the integrand for the one-loop matching 
factor, that is, (f(q)) = Z mj i(/ta). We find 



(/(q)log(qV)) 



2 5 
— log 2 (/ia) - — log(//.a) 

7T 37T 

1 5tt 

^ + T2-°- 821 



(33) 



where the numerical constant is the result of the numer- 
ical evaluation of the lattice moment (logarithms of mo 
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due to the anomalous dimension of the pole mass cancel 
in the scale for the matching coefficient). 

When \Z m .i([ia)\ is anomalously small, a proper eval- 
uation of q* requires the second-order expression 

log(«&) = «log(g 2 ))) ± [«log( 9 2 ))) 2 - «logV)))] 1/2 , 

(34) 

where the appropriate root (if the result is real) is usually 
made obvious by requiring continuity, and a physically 
reasonable value, for the resulting q* , as a function of 
the underlying parameters. The second moment of Z m ^\ 
is given by 

o in 

(/(g)logVa 2 )} = --log 3 (H log 2 (H 

(tt 2 -!}) 5 10 Ca 7tt 

+ — 5 log(Ma) - q o + 3 - 057 - ( 35 ) 

Sir 6ir 6tt 6b 

Results for the first- and (where appropriate) second- 
order scales, as functions of /ia, are shown in Fig. [S] A 
typical value is aq* — 1.877 at \ia = 1. 

III. RESULTS FOR LIGHT-QUARK MASSES 

The bare lattice masses for the strange and up/down 
quarks, on the MILC "coarse" and "fine" lattices, are 
given in Ref . |l3| . For the strange quark these are amo s — 
0.0390(l)(20j/u O c, and am 0s = 0.0272(1) (12)/u of , on 
the coarse and fine lattices respectively, where uq c — 
0.86774 and Uof — 0.87821 are tadpole normalization 
factors. The uncertainties are lattice statistical and sys- 
tematic errors, respectively, the latter due mainly to 
chiral extrapolation/interpolation. The lattice spacings 
can be found in Ref. 0, a^irse = 1-596(30) GeV, and 
a~^ e = 2.258(32) GeV. 

Following conventional practice we quote the light- 
quark MS masses at the scale /i = 2 GeV, taking three 
active flavors of quarks (rif = 3). The BLM scales on 
the two MILC lattices are then aq* oarso = 2.144 and 
aq^ ne = 1.752 j22|. This results in two-loop coefficients 
in Eq. of Z m , 2 | coarsc = 1.939(20), and Z m . 2 \ Rnc = 
2.270(20). We also require the couplings at the relevant 
scales, and for this purpose we use the recently deter- 
mined value a^ /=3) (7.5 GeV) = 0.2082(40) Q. We find 
H^co^e) = 0.2925(92), and a v (q* Rnc ) - 0.2713(76) 

Putting all of this together, we obtain the following 
values for the MS strange-quark mass, on the two lattices 

771^(2 GeV) = 84(5) MeV [coarse], 

771^(2 GeV) = 86(4) MeV [fine], (36) 

where the errors above are just from the simulation 
systematics; we note that these errors are correlated, 
because the two bare masses are obtained from a 
simultaneous-chiral fit to the two lattice spacings, which 



describes the effects of taste-changing and other dis- 
cretization effects in the staggered action ||. 

We consider continuum extrapolations of these values, 
based on the form of the expected leading-discretization 
errors, which are of 0(aya 2 ) (we find essentially identical 
results by assuming 0(a v a 2 ) errors 13]). In addition, 
we estimate the third-order perturbative correction to 
the matching factor, that is, we add a term Z m ^(fia)a v 
to the right-hand side of Eq. H26[) . and attempt to esti- 
mate its coefficient. To this end, we extend the lattice 
renormalization factor Eq. (JJJ to third order, including 
the logarithms from the three-loop anomalous dimension, 
which are fixed from the known MS expansion, along 
with the known third-order term in Eq. (J5J) for the MS 
mass j^, |2^. This leaves one unknown constant, in 
the lattice renormalization factor, A 30 in the notation 
of Eq. |0), which with the third-order logarithms deter- 
mines Z m ^(fia). In principle, one can extract A^q, and 
hence the third-order correction, from a simultaneous fit 
using bare lattice masses at several lattice spacings; this 
extends a technique first laid out in Ref. 9] (see also 
Ref. [13). 

With only the two available lattice spacings, and hav- 
ing also to include a discretization correction in the fit, 
we can only roughly bound the size of the next order in 
the perturbative expansion. We used constrained curve- 
fitting 0, 1.3 ll ] to include our expectation that the expan- 
sion is convergent (i.e., \Z m ^(fia)\ — O(l) in the notation 
of Eq. 

We tested this procedure by considering a fit to the 
second-order perturbative correction Z m ^,{\id), without 
a priori knowledge of the associated constant A20 : with 
the two lattice spacings as input, and including a dis- 
cretization correction, the fit returns A20 — 5.9 ± 1.9, in 
good agreement with Eq. (|25|l . The second-order fit also 
returns mf s (2 GeV) = (85 ± 11) MeV, in good agree- 
ment with the "bona-fide" two-loop values in Eq. 136|) 
[this also represents somewhat of an improvement com- 
pared to our earlier result |13| . which used only a pri- 
ori first-order perturbation theory, without a fit to the 
second-order correction, and which somewhat underesti- 
mated both the central value and the systematic error 
from the truncation of the perturbative series] . 

When this procedure is applied at third-order, with 
A 2 q input from Eq. (|2l>|l , the fit provides a reasonable es- 
timate of the relative systematic error on the MS mass, 
due to the third-order perturbative correction, of approx- 
imately 2 x a.y, or about 4% (there is no appreciable 
third-order correction to the MS mass, within this error). 
We also use these fits (which include a discretization cor- 
rection) to extract the central value of the strange-quark 
mass. Our final value is then 

m^(2 GeV) = 87(0) (4) (4) (0) MeV (37) 

where, following Ref. l ne respective errors are sta- 
tistical, lattice systematic, perturbative, and electromag- 
netic/isospin effects. 
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Our result for the ratio of the strange-quark mass to 
the up/down-quark masses is unchanged from Ref. [l3j . 
since the renormalization factor is mass independent, as 
we have verified explicitly in Sect. ITT1 through two-loops 
(and up to a negligible mass-dependent discretization 
correction) 

^=27.4(1)(4)(0)(1), (38) 

777 

where tti = \(m u + m^). Equivalently we have 

777^(2 GeV) = 3.2(0) (2) (2) (0) MeV (39) 

Using a recent determination of the ratio m u /rrid = 

0.43(0) (1)(0) (8) due to the MILC collaboration |H2, 
these results imply 

mf S (2 GeV) = 1.9(0) (1)(1) (2) MeV, 

777^(2 GeV) = 4.4(0) (2) (2) (2) MeV. (40) 

IV. DISCUSSION AND CONCLUSIONS 

Perturbation theory has once again shown itself to be 
an essential tool in high-precision phenomenological cal- 
culations from the lattice. The two-loop lattice diagrams 
for the renormalized-quark mass were conquered with 
a combination of algebraic and numerical techniques in 
this first-ever two-loop evaluation of a multiplicative "ki- 
netic" mass on the lattice. When combined with the 
known continuum matching from the pole mass to the 
MS mass, a very accurate determination of the light- 
quark masses was possible. The results presented here 
have a number of distinguishing features: two-loop per- 
turbation theory; rif = 2 + 1 simulations with two de- 
generate light quarks and a heavier strange quark; very 
small light quark masses from m s /8 to m s /2 which en- 
abled a partially quenched chiral fit with many terms 
to thousands of configurations; and extremely accurate 
determinations of the lattice spacings, which are equal 
within small errors when set from any of a wide variety 
of hadronic inputs. 

Most notable amongst our results is our new value for 
strange quark mass, mf (2 GeV) = 87(0) (4) (4) (0) MeV, 
where the respective errors are lattice statistical, lat- 
tice systematic (mostly due to the chiral extrapola- 
tion/interpolation), perturbative, and due to electromag- 
netic/isospin effects. The two-loop matching has in- 
creased the central value of our estimates of the light- 
quark masses with respect to our previous one-loop de- 
termination 0| by about 1.5 standard deviations, based 
on the previous estimate of the perturbation-theory un- 
certainty. The systematic uncertainty from perturbation 
theory has been reduced by about a factor of two, and is 
now only about 4%, the same size as the current lattice 
systematic uncertainty, the latter due mainly to the chi- 
ral extrapolation/interpolation. We anticipate that the 



present estimate of the perturbative uncertainty could 
be reduced somewhat further, if additional lattice spac- 
ings become available, by using the NNLO perturbation 
theory presented here to improve the estimate of the 
third-order perturbative correction, along the lines that 
we have implemented above. 

The strange-quark mass determination has histori- 
cally generated some controversy, with widely different 
values having been obtained from different approaches. 
This is reflected in the large uncertainty in the Particle 
Data Group's most recent best estimate, (2 GeV) = 
(105 ± 25) MeV js3l; our result represents a significant 
improvement in precision, resulting from an aggressive 
effort to understand and reduce all sources of systematic 
error. 

An obvious advantage of our result is that it has been 
obtained with the correct description of the sea, that is, 
with rif = 2 + 1 flavors of dynamical quarks. There is 
only one other three-flavor result, which is due to the 
CP-PACS and JLQCD collaborations (which did simula- 
tions at much larger u/d quark masses than in the MILC 
"asqtad" ensembles), which recently reported a value of 

mp(2 GeV) = 87(6) MeV [H; however they do not 
include a full error analysis, and in particular the error 
from missing higher orders in the perturbative matching 
alone should be comparable to that in our older result, 
and significantly larger than the errors that we report 
here with a 2-loop analysis. 

It appears that the most recent estimates of the 
strange-quark mass extracted from simulations with only 
two flavors of sea quarks are systematically higher than 
the estimates with the correct 77/, although the other 
systematic errors are too large to allow for a definitive 
assessment (noting that these two-flavor determinations 
were also done with different definitions of the quark 
mass and determinations of the lattice spacing from 7*0 
using different physical values for 7"o). The two-flavor de- 
termination from the QCDSF-UKQCD collaboration is 
777^(2 GeV) = 100 - 130 MeV H3, the ALPHA collab- 
oration value is 97(22) MeV [3(|, and the Rome value is 
101(8) (^lg) MeV [33. In this connection, we analyzed 
quenched simulations of the "asqtad" action by the MILC 
collaboration |3J, |5j , and we find that this also leads to 
a somewhat larger value of the strange quark mass, of 
about 96 MeV quenched (using simple linear interpola- 
tions of the MILC quenched-meson spectrum results, and 
our two-loop perturbative matching formula at 77/ = 0). 

We are currently in the process of applying our NNLO 
matching calculation to heavy-quark masses, in order to 
complete a high-precision determination of the funda- 
mental parameters of QCD. Work is also underway on 
the NNLO matching calculations for important hadronic 
matrix elements, especially those such as the decay con- 
stants fn and fs and other form-factors of particular 
relevance to heavy-flavor physics. 



10 



Acknowledgments 

This work was supported by the US Department of 
Energy, the US National Science Foundation, the Natural 



Science and Engineering Research Council of Canada, the 
Particle Physics and Astronomy Research Council of the 
UK, and Barclays Capital. We thank Matthew Nobes 
and Kent Hornbostel for fruitful discussions. 



[1] G. P. Lepage and P. B. Mackenzie, Phys. Rev. D 48, 2250 
(1993) 

[2] G. P. Lepage, Phys. Rev. D 59, 074502 (1999). 

[3] K. Orginos and D. Toussaint, Phys. Rev. D 59, 014501 

(1999); C. W. Bernard et al, MILC Collaboration, Phys. 

Rev. D 64, 054506 (2001) 
[4] W. J. Lee and S. R. Sharpe, Phys. Rev. D 60, 114503 

(1999); C. Aubin and C. Bernard, Phys. Rev. D 68, 

034014 (2003) 

[5] C. Aubin et al, MILC Collaboration, Phys. Rev. D 70, 
094505 (2004). 

[6] C. T. H. Davies et al, HPQCD, Fermilab and MILC 
Collaborations, Phys. Rev. Lett. 92 022001 (2004). 

[7] H. D. Trottier, Nucl. Phys. Proc. Suppl. 129, 142 (2004). 

[8] Q. J. Mason, Ph.D. thesis, Cornell University (2004). 

[9] Q. Mason et al, HPQCD Collaboration, Phys. Rev. Lett. 
95 (2005) 052002. 
[10] E. Follana and H. Panagopoulos, Phys. Rev. D 63 (2000) 
017501. 

[11] U. M. Heller and F. Karsch, Nucl. Phys. B 251, 254 
(1985); G. Martinelli and C. T. Sachrajda, Nucl. Phys. 
B 559, 429 (1999). 

[12] A. Gray et al, HPQCD Collaboration, |hep-lat /0507013| 

[13] C. Aubin et al, HPQCD, MILC, and UKQCD Collabo- 
rations, Phys. Rev. D. 70 (2004) 031504. 

[14] J. Hein, Q. Mason, G. P. Lepage and H. Trottier, Nucl. 
Phys. Proc. Suppl. 106, 236 (2002). 

[15] W. J. Lee and S. R. Sharpe, Phys. Rev. D 66, 114501 
(2002). 

[16] Q. Mason and H. D. Trottier, in preparation. 

[17] E. Follana, A. Hart and C. T. H. Davies, Phys. Rev. Lett. 
93, 241601 (2004); K. Y. Wong and R. M. Woloshyn, 
Phys. Rev. D 71, 094508 (2005); S. Durr, C. Hoelbling 
and U. Wenger, Phys. Rev. D 70, 094502 (2004). 

[18] Q. Mason et al, HPQCD collaboration, Nucl. Phys. 
Proc. Suppl. 119, 446 (2003); E. Follana et al, HPQCD 
Collaboration, ibid., 129 & 130, 384 (2004). 

[19] R. Tarrach, Nucl. Phys. B 183, 384 (1981). 

[20] N. Gray, D.J. Broadhurst, W. Grafe and K. Schilcher, Z. 



Phys. C 48, 673 (1990); D. J. Broadhurst, N. Gray and 
K. Schilcher, ibid. 52, 111 (1991). 
[21] K. G. Chetyrkin and M. Steinhauser, Nucl. Phys. B 573, 
617 (2000). 

[22] K. Melnikov and T. v. Ritbergen, Phys. Lett. B 482 
(2000) 99. 

[23] Q. Mason, H. D. Trottier, R. Horgan, hep-lat/0510053 
[24] Y. Schroder, Phys. Lett. B 447, 321 (1999) 
[25] A. S. Kronfeld, Phys. Rev. D 58, 051501 (1998). 
[26] G. P. Lepage, J. Comput. Phys. 27 (1978) 192. 
[27] For a review see, e.g., M. Beneke, Phys. Rept. 317, 1 
(1999). 

[28] S. J. Brodsky, G. P. Lepage and P. B. Mackenzie, Phys. 
Rev. D 28, 228 (1983); K. Hornbostel, G. P. Lepage and 
C. Morningstar, Phys. Rev. D 67, 034023 (2003). 

[29] The values of the q* scales used here are more accu- 
rate than the slightly larger values given in Ref. |l3j|. 
where an approximate determination was made of the 
average-loop momentum circulating in the continuum 
self-energy. In addition an earlier, approximate (and 

somewhat smaller), determination of Oy f 3 '(7.5 GeV) 

was employed in Ref. [13J. 
[30] H. D. Trottier, N. H. Shakespeare, G. P. Lepage and 

P. B. Mackenzie, Phys. Rev. D 65, 094502 (2002). 
[31] G. P. Lepage et al, Nucl. Phys. Proc. Suppl. 106, 12 

(2002) 

[32] C. Aubin et al, MILC collaboration, Phys. Ref. D 70, 
114501 (2004). 

[33] S. Eidelman et al, Particle Data Group, Phys. Lett. B 
592, 1 (2004). 

[34] T. Ishikawa et al , CP-PACS and JLQCD collaborations, 

hep-lat/0509142 
[35j M. Gockeler e t al, QCDSF-UKQCD collaboration, 

hep-lat/0509159 
[36] M. Delia Mort e et al, ALPHA collaboration, 

|hep-lat/0509073"l 
[37] D. Bccircvic et al, hep-lat/0509091 



